n.cntrys = length(unique(trustleg$Country))
n.yrs = 2020-year0leg # estimates up to 2020
n.resp = dim(trustleg)[1]
n.itm.cnt = length(unique(trustleg$ItemCnt))
cntrys = as.numeric(factor(trustleg$Country))
cnt.names = levels(trustleg$Country)
items = as.numeric(factor(trustleg$Item))
yrs = trustleg$Year
itm.cnts = as.numeric(factor(trustleg$ItemCnt))
mean.resp.log = logit(mean(trustleg$Response))
# create item-country length indicator for items
item.ind.kp = rep(0, length(levels(trustleg$ItemCnt)))
for(i in 1:length(levels(trustleg$Item))) {
item.ind.kp[grepl(levels(trustleg$Item)[i], levels(trustleg$ItemCnt))] = i
}
item.ind.len = sapply(lapply(levels(trustleg$Item), function(x) grep(x, levels(trustleg$ItemCnt))), length)
## Fit stan model
# specify data for stan
dat.1 = list(N=n.resp, K=n.items, R=n.yrs, J=n.cntrys, P=n.itm.cnt, jj=cntrys, rr=yrs,
pp=itm.cnts, kk=items, it_len=item.ind.len,
x=trustleg$RespN, samp=trustleg$Sample, mn_resp_log=mean.resp.log)
sapply(dat.1, summary)
# parameters to save
pars.1 = c("Sigma","Omega","sigma_delta","sigma_theta","phi","mu_lambda","lambda","gamm","delta",
"theta","x_pred","log_lik")
# iterations for MCMC simulation
n.iter = 1000
n.warm = 500
n.samp = n.iter - n.warm
n.chn = 4
# compile model
stan.mod = cmdstan_model('stan_mod6_v6_rep.stan')
# compile model
stan.mod = cmdstan_model('stan.mod6.7_newstanrcode.stan')
# Stan fit
fit.mod= stan.mod$sample(
data = dat.1,
chains = n.chn,
init = 0.1,
parallel_chains = n.chn,
iter_warmup = n.warm,
iter_sampling = n.samp,
refresh = round(n.iter/20, 0),
adapt_delta = 0.8,
max_treedepth = 15,
save_warmup = FALSE
)
## Check convergence
# Examine model fit
res = fit.mod
res$cmdstan_diagnose()
res.tab = res$print(pars.1, max_rows=80, digits=3)
sum = res$summary(pars.1)
print(sum[order(sum$rhat, decreasing=TRUE), ], n=50)
res_neff_ratio = neff_ratio(res)
res_neff_ratio[order(res_neff_ratio, decreasing=FALSE)][1:50]
# traceplot
tp.pars = c("Sigma[1,1]","Sigma[2,2]","Omega[1,2]","sigma_theta","sigma_delta","mu_lambda",
"phi","delta[23]","theta[31,121]")
tp = bayesplot::mcmc_trace(res$draws(tp.pars), size=0.3, np=nuts_params(res))
tp
## Extract and save mood estimates
theta.m.out = apply(res$draws("theta"), 3, as.vector)
(theta.m.mean = mean(as.vector(theta.m.out)))
(theta.m.sd = sd(as.vector(theta.m.out)))
theta.m.std = (theta.m.out - theta.m.mean) / theta.m.sd # standardize
theta.m.t = apply(theta.m.std, 1, function(x) t(x) )
theta.pe = apply(theta.m.t, 1, mean)
theta.u95 = apply(theta.m.t, 1, quantile, probs=c(0.975))
theta.l95 = apply(theta.m.t, 1, quantile, probs=c(0.025))
theta.sd = apply(theta.m.t, 1, sd)
theta.m.df = data.frame(Country=rep(cnt.names, each=n.yrs),
Year=rep(1973:2020, times=n.cntrys), leg=theta.pe,
leg_u95=theta.u95, leg_l95=theta.l95, leg_sd=theta.sd)
# remove estimates before first survey year and create a trimmed dataset
first.yr = data.frame(Country = levels(trustleg$Country),
First_yr = as.vector(
by(trustleg, trustleg$Country, function(x) min(as.numeric(x$Year)) + year0leg)))
theta.trim = merge(theta.m.df, first.yr, by="Country", all.x=TRUE)
cnts = theta.trim[theta.trim$Year==2008, "Country"]
frst.yr = theta.trim[theta.trim$Year==2008, "First_yr"]
theta.trim$leg_trim = theta.trim$leg
theta.trim$leg_trim = ifelse(theta.trim$Year < theta.trim$First_yr, NA, theta.trim$leg_trim)
theta.trim = theta.trim[order(theta.trim$Country, theta.trim$Year), ]
theta.trim = theta.trim[!is.na(theta.trim$leg_trim),]
theta.trim$leg_trim = NULL
View(theta.trim)
# save country-year point estimates
write.csv(theta.trim, "leg_mood_est.csv", row.names=FALSE)
###### Estimating Trust in Parliament, 1973-2022
###### Running Stan Models Using cmdstanr
###### v5
# libraries
library(arm)
library(dplyr)
library(tidyr)
library(loo)
library(ggplot2)
library(bayesplot)
library(cmdstanr)
library(posterior)
library(rio)
library(tidyverse)
setwd("C:/Users/ba72loko/projects/Democratic Resilience (with Aurel)/data/political trust BJPS study")
# options
options("cmdstanr_verbose" = TRUE)
options(mc.cores = parallel::detectCores())
# read trust data
trustparl = read.csv("trends_parl.csv")
## Edit data
# remove NAs
trustparl = trustparl[!trustparl$Response==0, ]
# order
trustparl = arrange(trustparl, Country, Year)
# set first year
year0parl = 1972 # = year before first year of available survey data
trustparl = trustparl[trustparl$Year > year0parl,]
# I originally used the item variable with wording variation as the item variable for this analysis
# This resulted in many item-country combinations only including one observation and thus there were very high rhat values
# Try with a broader item variable instead, which only uses the survey source (as Chris Claassen does)
trustparl <- dplyr::select(trustparl, -Item)
trustparl <- dplyr::rename(trustparl, Item = Project)
# Recode the value of Item to make the EB and AsiaBarometer data source identifiers more specific
# Initially, the grep code for item-country length identifiers below thought that many other observations were also from these sources
trustparl$Item[trustparl$Item=="EB"] <- "SSEB" # stands for "Standard and Special Eurobarometer
trustparl$Item[trustparl$Item=="AsB"] <- "AsiaB"
# create item by country indicators
trustparl = unite(trustparl, ItemCnt, c(Item, Country), sep = "_", remove = FALSE)
# create item by region indicators
trustparl = unite(trustparl, ItemReg, c(Item, regpol6), sep = "_", remove = FALSE)
# identify countries with few years of data
cnt.obs.years = rowSums(table(trustparl$Country, trustparl$Year) > 0)
sort(cnt.obs.years)
# run the next line to drop countries with less than 2 years of data
trustparl = trustparl[trustparl$Country %in% levels(factor(trustparl$Country))[cnt.obs.years > 1], ]
length(unique(trustparl$Country))
## Prepare data for stan
# factorise
trustparl$Country = as.factor(as.character(trustparl$Country))
trustparl$Item = as.factor(as.character(trustparl$Item))
trustparl$ItemCnt = as.factor(as.character(trustparl$ItemCnt))
trustparl$Year = trustparl$Year-year0parl
# extract data
n.items = length(unique(trustparl$Item))
n.cntrys = length(unique(trustparl$Country))
n.yrs = 2020-year0parl # estimates up to 2020
n.resp = dim(trustparl)[1]
n.itm.cnt = length(unique(trustparl$ItemCnt))
cntrys = as.numeric(factor(trustparl$Country))
cnt.names = levels(trustparl$Country)
items = as.numeric(factor(trustparl$Item))
yrs = trustparl$Year
itm.cnts = as.numeric(factor(trustparl$ItemCnt))
mean.resp.log = logit(mean(trustparl$Response))
# create item-country length indicator for items
item.ind.kp = rep(0, length(levels(trustparl$ItemCnt)))
for(i in 1:length(levels(trustparl$Item))) {
item.ind.kp[grepl(levels(trustparl$Item)[i], levels(trustparl$ItemCnt))] = i
}
item.ind.len = sapply(lapply(levels(trustparl$Item), function(x) grep(x, levels(trustparl$ItemCnt))), length)
## Fit stan model
# specify data for stan
dat.1 = list(N=n.resp, K=n.items, R=n.yrs, J=n.cntrys, P=n.itm.cnt, jj=cntrys, rr=yrs,
pp=itm.cnts, kk=items, it_len=item.ind.len,
x=trustparl$RespN, samp=trustparl$Sample, mn_resp_log=mean.resp.log)
sapply(dat.1, summary)
# parameters to save
pars.1 = c("Sigma","Omega","sigma_delta","sigma_theta","phi","mu_lambda","lambda","gamm","delta",
"theta","x_pred","log_lik")
# iterations for MCMC simulation
n.iter = 1000
n.warm = 500
n.samp = n.iter - n.warm
n.chn = 4
# compile model
stan.mod = cmdstan_model('stan.mod6.7_newstanrcode.stan')
# Stan fit
fit.mod= stan.mod$sample(
data = dat.1,
chains = n.chn,
init = 0.1,
parallel_chains = n.chn,
iter_warmup = n.warm,
iter_sampling = n.samp,
refresh = round(n.iter/20, 0),
adapt_delta = 0.8,
max_treedepth = 15,
save_warmup = FALSE
)
## Check convergence
# Examine model fit
res = fit.mod
res$cmdstan_diagnose()
res.tab = res$print(pars.1, max_rows=80, digits=3)
sum = res$summary(pars.1)
print(sum[order(sum$rhat, decreasing=TRUE), ], n=50)
res_neff_ratio = neff_ratio(res)
res_neff_ratio[order(res_neff_ratio, decreasing=FALSE)][1:50]
# traceplot
tp.pars = c("Sigma[1,1]","Sigma[2,2]","Omega[1,2]","sigma_theta","sigma_delta","mu_lambda",
"phi","delta[23]","theta[31,121]")
tp = bayesplot::mcmc_trace(res$draws(tp.pars), size=0.3, np=nuts_params(res))
tp
## Extract and save mood estimates
theta.m.out = apply(res$draws("theta"), 3, as.vector)
(theta.m.mean = mean(as.vector(theta.m.out)))
(theta.m.sd = sd(as.vector(theta.m.out)))
theta.m.std = (theta.m.out - theta.m.mean) / theta.m.sd # standardize
theta.m.t = apply(theta.m.std, 1, function(x) t(x) )
theta.pe = apply(theta.m.t, 1, mean)
theta.u95 = apply(theta.m.t, 1, quantile, probs=c(0.975))
theta.l95 = apply(theta.m.t, 1, quantile, probs=c(0.025))
theta.sd = apply(theta.m.t, 1, sd)
theta.m.df = data.frame(Country=rep(cnt.names, each=n.yrs),
Year=rep(1973:2020, times=n.cntrys), parl=theta.pe,
parl_u95=theta.u95, parl_l95=theta.l95, parl_sd=theta.sd)
# remove estimates before first survey year and create a trimmed dataset
first.yr = data.frame(Country = levels(trustparl$Country),
First_yr = as.vector(
by(trustparl, trustparl$Country, function(x) min(as.numeric(x$Year)) + year0parl)))
theta.trim = merge(theta.m.df, first.yr, by="Country", all.x=TRUE)
cnts = theta.trim[theta.trim$Year==2008, "Country"]
frst.yr = theta.trim[theta.trim$Year==2008, "First_yr"]
theta.trim$parl_trim = theta.trim$parl
theta.trim$parl_trim = ifelse(theta.trim$Year < theta.trim$First_yr, NA, theta.trim$parl_trim)
theta.trim = theta.trim[order(theta.trim$Country, theta.trim$Year), ]
theta.trim = theta.trim[!is.na(theta.trim$parl_trim),]
theta.trim$parl_trim = NULL
View(theta.trim)
# save country-year point estimates
write.csv(theta.trim, "parl_mood_est.csv", row.names=FALSE)
###### Estimating Trust in the police
###### Running Stan Models Using cmdstanr
###### v5
# libraries
library(arm)
library(dplyr)
library(tidyr)
library(loo)
library(ggplot2)
library(bayesplot)
library(cmdstanr)
library(posterior)
library(rio)
library(tidyverse)
setwd("C:/Users/ba72loko/projects/Democratic Resilience (with Aurel)/data/political trust BJPS study")
# options
options("cmdstanr_verbose" = TRUE)
options(mc.cores = parallel::detectCores())
# read trust data
trustpolice = read.csv("trends_police.csv")
## Edit data
# remove NAs
trustpolice = trustpolice[!trustpolice$Response==0, ]
# order
trustpolice = arrange(trustpolice, Country, Year)
# set first year
year0police = 1980 # = year before first year of available survey data
trustpolice = trustpolice[trustpolice$Year > year0police,]
# I originally used the item variable with wording variation as the item variable for this analysis
# This resulted in many item-country combinations only including one observation and thus there were very high rhat values
# Try with a broader item variable instead, which only uses the survey source (as Chris Claassen does)
trustpolice <- dplyr::select(trustpolice, -Item)
trustpolice <- dplyr::rename(trustpolice, Item = Project)
# Recode the value of Item to make the EB and AsiaBarometer data source identifiers more specific
# Initially, the grep code for item-country length identifiers below thought that many other observations were also from these sources
trustpolice$Item[trustpolice$Item=="EB"] <- "SSEB" # stands for "Standard and Special Eurobarometer
trustpolice$Item[trustpolice$Item=="AsB"] <- "AsiaB" # stands for "Standard and Special Eurobarometer
# create item by country indicators
trustpolice = unite(trustpolice, ItemCnt, c(Item, Country), sep = "_", remove = FALSE)
# create item by region indicators
trustpolice = unite(trustpolice, ItemReg, c(Item, regpol6), sep = "_", remove = FALSE)
# identify countries with few years of data
cnt.obs.years = rowSums(table(trustpolice$Country, trustpolice$Year) > 0)
sort(cnt.obs.years)
# run the next line to drop countries with less than 2 years of data
trustpolice = trustpolice[trustpolice$Country %in% levels(factor(trustpolice$Country))[cnt.obs.years > 1], ]
length(unique(trustpolice$Country))
## Prepare data for stan
# factorise
trustpolice$Country = as.factor(as.character(trustpolice$Country))
trustpolice$Item = as.factor(as.character(trustpolice$Item))
trustpolice$ItemCnt = as.factor(as.character(trustpolice$ItemCnt))
trustpolice$Year = trustpolice$Year-year0police
# extract data
n.items = length(unique(trustpolice$Item))
n.cntrys = length(unique(trustpolice$Country))
n.yrs = 2020-year0police # estimates up to 2020
n.resp = dim(trustpolice)[1]
n.itm.cnt = length(unique(trustpolice$ItemCnt))
cntrys = as.numeric(factor(trustpolice$Country))
cnt.names = levels(trustpolice$Country)
items = as.numeric(factor(trustpolice$Item))
yrs = trustpolice$Year
itm.cnts = as.numeric(factor(trustpolice$ItemCnt))
mean.resp.log = logit(mean(trustpolice$Response))
# create item-country length indicator for items
item.ind.kp = rep(0, length(levels(trustpolice$ItemCnt)))
for(i in 1:length(levels(trustpolice$Item))) {
item.ind.kp[grepl(levels(trustpolice$Item)[i], levels(trustpolice$ItemCnt))] = i
}
item.ind.len = sapply(lapply(levels(trustpolice$Item), function(x) grep(x, levels(trustpolice$ItemCnt))), length)
## Fit stan model
# specify data for stan
dat.1 = list(N=n.resp, K=n.items, R=n.yrs, J=n.cntrys, P=n.itm.cnt, jj=cntrys, rr=yrs,
pp=itm.cnts, kk=items, it_len=item.ind.len,
x=trustpolice$RespN, samp=trustpolice$Sample, mn_resp_log=mean.resp.log)
sapply(dat.1, summary)
# parameters to save
pars.1 = c("Sigma","Omega","sigma_delta","sigma_theta","phi","mu_lambda","lambda","gamm","delta",
"theta","x_pred","log_lik")
# iterations for MCMC simulation
n.iter = 1000
n.warm = 500
n.samp = n.iter - n.warm
n.chn = 4
# compile model
stan.mod = cmdstan_model('stan.mod6.7_newstanrcode.stan')
# Stan fit
fit.mod= stan.mod$sample(
data = dat.1,
chains = n.chn,
init = 0.1,
parallel_chains = n.chn,
iter_warmup = n.warm,
iter_sampling = n.samp,
refresh = round(n.iter/20, 0),
adapt_delta = 0.8,
max_treedepth = 15,
save_warmup = FALSE
)
## Check convergence
# Examine model fit
res = fit.mod
res$cmdstan_diagnose()
res.tab = res$print(pars.1, max_rows=80, digits=3)
sum = res$summary(pars.1)
print(sum[order(sum$rhat, decreasing=TRUE), ], n=50)
res_neff_ratio = neff_ratio(res)
res_neff_ratio[order(res_neff_ratio, decreasing=FALSE)][1:50]
# traceplot
tp.pars = c("Sigma[1,1]","Sigma[2,2]","Omega[1,2]","sigma_theta","sigma_delta","mu_lambda",
"phi","delta[23]","theta[31,121]")
tp = bayesplot::mcmc_trace(res$draws(tp.pars), size=0.3, np=nuts_params(res))
tp
## Extract and save mood estimates
theta.m.out = apply(res$draws("theta"), 3, as.vector)
(theta.m.mean = mean(as.vector(theta.m.out)))
(theta.m.sd = sd(as.vector(theta.m.out)))
theta.m.std = (theta.m.out - theta.m.mean) / theta.m.sd # standardize
theta.m.t = apply(theta.m.std, 1, function(x) t(x) )
theta.pe = apply(theta.m.t, 1, mean)
theta.u95 = apply(theta.m.t, 1, quantile, probs=c(0.975))
theta.l95 = apply(theta.m.t, 1, quantile, probs=c(0.025))
theta.sd = apply(theta.m.t, 1, sd)
theta.m.df = data.frame(Country=rep(cnt.names, each=n.yrs),
Year=rep(1981:2020, times=n.cntrys), police=theta.pe,
police_u95=theta.u95, police_l95=theta.l95, police_sd=theta.sd)
# remove estimates before first survey year and create a trimmed dataset
first.yr = data.frame(Country = levels(trustpolice$Country),
First_yr = as.vector(
by(trustpolice, trustpolice$Country, function(x) min(as.numeric(x$Year)) + year0police)))
theta.trim = merge(theta.m.df, first.yr, by="Country", all.x=TRUE)
cnts = theta.trim[theta.trim$Year==2008, "Country"]
frst.yr = theta.trim[theta.trim$Year==2008, "First_yr"]
theta.trim$police_trim = theta.trim$police
theta.trim$police_trim = ifelse(theta.trim$Year < theta.trim$First_yr, NA, theta.trim$police_trim)
theta.trim = theta.trim[order(theta.trim$Country, theta.trim$Year), ]
theta.trim = theta.trim[!is.na(theta.trim$police_trim),]
theta.trim$police_trim = NULL
# save country-year point estimates
write.csv(theta.trim, "police_mood_est.csv", row.names=FALSE)
###### Estimating Trust in Political Parties, 1990-2022
###### Running Stan Models Using cmdstanr
###### v5
# libraries
library(arm)
library(dplyr)
library(tidyr)
library(loo)
library(ggplot2)
library(bayesplot)
library(cmdstanr)
library(purrr)
library(rio)
library(tidyverse)
setwd("C:/Users/ba72loko/projects/Democratic Resilience (with Aurel)/data/political trust BJPS study")
# options
options("cmdstanr_verbose" = TRUE)
options(mc.cores = parallel::detectCores())
# read trust data
trustpolpar = read.csv("trends_polpar.csv")
## Edit data
# remove NAs
trustpolpar = trustpolpar[!trustpolpar$Response==0, ]
# order
trustpolpar = arrange(trustpolpar, Country, Year)
# I originally used the item variable with wording variation as the item variable for this analysis
# Use a broader item variable instead, which only uses the survey source (as Chris Claassen does), generally results in more dynamic (less flat) trends
trustpolpar <- dplyr::select(trustpolpar, -Item)
trustpolpar <- dplyr::rename(trustpolpar, Item = Project)
# Recode the value of Item to make the EB and AsiaBarometer data source identifiers more specific
# Initially, the grep code for item-country length identifiers below thought that many other observations were also from these sources
trustpolpar$Item[trustpolpar$Item=="EB"] <- "SSEB" # stands for "Standard and Special Eurobarometer
trustpolpar$Item[trustpolpar$Item=="AsB"] <- "AsiaB" # stands for "Standard and Special Eurobarometer
# set first year
year0polpar = 1989 # = year before first year of available survey data
trustpolpar = trustpolpar[trustpolpar$Year > year0polpar,]
# create item by country indicators
trustpolpar = unite(trustpolpar, ItemCnt, c(Item, Country), sep = "_", remove = FALSE)
# create item by region indicators
trustpolpar = unite(trustpolpar, ItemReg, c(Item, regpol6), sep = "_", remove = FALSE)
# identify countries with few years of data
cnt.obs.years = rowSums(table(trustpolpar$Country, trustpolpar$Year) > 0)
sort(cnt.obs.years)
# run the next line to drop countries with less than 2 years of data
trustpolpar = trustpolpar[trustpolpar$Country %in% levels(factor(trustpolpar$Country))[cnt.obs.years > 1], ]
length(unique(trustpolpar$Country))
## Prepare data for stan
# factorise
trustpolpar$Country = as.factor(as.character(trustpolpar$Country))
trustpolpar$Item = as.factor(as.character(trustpolpar$Item))
trustpolpar$ItemCnt = as.factor(as.character(trustpolpar$ItemCnt))
trustpolpar$Year = trustpolpar$Year-year0polpar
# extract data
n.items = length(unique(trustpolpar$Item))
n.cntrys = length(unique(trustpolpar$Country))
n.yrs = 2020-year0polpar # estimates up to 2020
n.resp = dim(trustpolpar)[1]
n.itm.cnt = length(unique(trustpolpar$ItemCnt))
cntrys = as.numeric(factor(trustpolpar$Country))
cnt.names = levels(trustpolpar$Country)
items = as.numeric(factor(trustpolpar$Item))
yrs = trustpolpar$Year
itm.cnts = as.numeric(factor(trustpolpar$ItemCnt))
mean.resp.log = logit(mean(trustpolpar$Response))
# create item-country length indicator for items
item.ind.kp = rep(0, length(levels(trustpolpar$ItemCnt)))
for(i in 1:length(levels(trustpolpar$Item))) {
item.ind.kp[grepl(levels(trustpolpar$Item)[i], levels(trustpolpar$ItemCnt))] = i
}
item.ind.len = sapply(lapply(levels(trustpolpar$Item), function(x) grep(x, levels(trustpolpar$ItemCnt))), length)
## Fit stan model
# specify data for stan
dat.1 = list(N=n.resp, K=n.items, R=n.yrs, J=n.cntrys, P=n.itm.cnt, jj=cntrys, rr=yrs,
pp=itm.cnts, kk=items, it_len=item.ind.len,
x=trustpolpar$RespN, samp=trustpolpar$Sample, mn_resp_log=mean.resp.log)
sapply(dat.1, summary)
# parameters to save
pars.1 = c("Sigma","Omega","sigma_delta","sigma_theta","phi","mu_lambda","lambda","gamm","delta",
"theta","x_pred","log_lik")
# iterations for MCMC simulation
n.iter = 1000
n.warm = 500
n.samp = n.iter - n.warm
n.chn = 4
# compile model
stan.mod = cmdstan_model('stan.mod6.7_newstanrcode.stan')
# Stan fit
fit.mod= stan.mod$sample(
data = dat.1,
chains = n.chn,
init = 0.1,
parallel_chains = n.chn,
iter_warmup = n.warm,
iter_sampling = n.samp,
refresh = round(n.iter/20, 0),
adapt_delta = 0.99,
max_treedepth = 13,
save_warmup = FALSE
)
## Check convergence
# Examine model fit
res = fit.mod
res$cmdstan_diagnose()
res.tab = res$print(pars.1, max_rows=80, digits=3)
sum = res$summary(pars.1)
print(sum[order(sum$rhat, decreasing=TRUE), ], n=50)
res_neff_ratio = neff_ratio(res)
res_neff_ratio[order(res_neff_ratio, decreasing=FALSE)][1:50]
# traceplot
tp.pars = c("Sigma[1,1]","Sigma[2,2]","Omega[1,2]","sigma_theta","sigma_delta","mu_lambda",
"phi","delta[23]")
tp = bayesplot::mcmc_trace(res$draws(tp.pars), size=0.3, np=nuts_params(res))
tp
## Extract and save mood estimates
theta.m.out = apply(res$draws("theta"), 3, as.vector)
(theta.m.mean = mean(as.vector(theta.m.out)))
(theta.m.sd = sd(as.vector(theta.m.out)))
theta.m.std = (theta.m.out - theta.m.mean) / theta.m.sd # standardize
theta.m.t = apply(theta.m.std, 1, function(x) t(x) )
theta.pe = apply(theta.m.t, 1, mean)
theta.u95 = apply(theta.m.t, 1, quantile, probs=c(0.975))
theta.l95 = apply(theta.m.t, 1, quantile, probs=c(0.025))
theta.sd = apply(theta.m.t, 1, sd)
theta.m.df = data.frame(Country=rep(cnt.names, each=n.yrs),
Year=rep(1990:2020, times=n.cntrys), polpar=theta.pe,
polpar_u95=theta.u95, polpar_l95=theta.l95, polpar_sd=theta.sd)
# remove estimates before first survey year and create a trimmed dataset
first.yr = data.frame(Country = levels(trustpolpar$Country),
First_yr = as.vector(
by(trustpolpar, trustpolpar$Country, function(x) min(as.numeric(x$Year)) + year0polpar)))
theta.trim = merge(theta.m.df, first.yr, by="Country", all.x=TRUE)
cnts = theta.trim[theta.trim$Year==2008, "Country"]
frst.yr = theta.trim[theta.trim$Year==2008, "First_yr"]
theta.trim$polpar_trim = theta.trim$polpar
theta.trim$polpar_trim = ifelse(theta.trim$Year < theta.trim$First_yr, NA, theta.trim$polpar_trim)
theta.trim = theta.trim[order(theta.trim$Country, theta.trim$Year), ]
theta.trim = theta.trim[!is.na(theta.trim$polpar_trim),]
theta.trim$polpar_trim = NULL
# save country-year point estimates
write.csv(theta.trim, "polpar_mood_est.csv", row.names=FALSE)
